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FOREWORD 
(Nontechnical  summary) 

An  improvement  in  NUMIT-ONE  (AFRRI  Scientific  Report  SR65-1,  October 
1965)  is  described.  Simpson's  Rule  is  still  used  to  perform  the  integration  but 
changes  in  the  coding  remove  the  need  for  a  large  array.  The  number  of  multiplica¬ 
tion  operations  done  during  a  run  has  also  been  vastly  reduced  resulting  in  significant 
savings  in  machine  time. 
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ABSTRACT 

An  improvement  on  NUMIT-ONE,  a  one-dimensional  numerical  integration  pro¬ 
gram  for  the  IBM  1620,  is  reported.  Accuracy,  input  and  output  formats,  source 
program,  flow  charts  and  test  problems  are  discussed. 
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I.  INTRODUCTION 


This  is  an  amendment  to  NUMIT-ONE,  a  numerical  integration  program  in  one- 
dimension  for  the  IBM  1620  computer,  AFRRI  Scientific  Report  SR65-1.  Appendix  D 
contains  a  glossary  of  terms  used  in  this  report. 

Further  analysis  of  the  problem  resulted  in  the  elimination  of  a  large  array 
thereby  making  the  program  smaller  and  the  number  of  multiplication  operations 
fewer.  This  reduction  in  size  enabled  compilation  with  a  larger  precision  parameter* 
(thereby  increasing  accuracy)  and  the  fewer  operations  decreased  the  length  of  running 
time. 


Successive  approximations  to  the  integral  are  now  computed  with  J  successively 


1,  2,  . . . ,  using 


AXi 


Fj  (a,b)  =  -3^  [f  (a)  +  f (b)  +  2  .  Cj  +  4  .  Dj] 


in  which 


AXj 


|  b-a  I 

I  2  J  I 


2J-l 
2  1 


Cj 


Dj 


£  f(a  +  2m  •  AXj),  J>1;  Cj  =  0  if  J  =  1 
m=l 


2J 

2 

£  f(a  +  [2m-l]  AXj) 
m=l 


Thus  for  the  J  +  1th  iteration: 


Cj+l  =  Dj+Cj 

*  A  precision  parameter  is  the  number  of  arithmetic  places  to  which  quantities  will  be 
held  during  computation  and  may  be  specified  on  a  control  card  that  precedes  the 
rest  of  the  deck  in  IBM  1620  FORTRAN  II. 
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and  only  two  summed  quantities  in  the  main  equation  need  be  multiplied  rather  than 


entire  arrays. 

Input  and  output  formats  are  slightly  changed  from  the  original  version. 

II.  ACCURACY 

Reduced  truncation  error  has  improved  the  accuracy.  An  inherent  weakness  in 
Simpson's  Rule  makes  it  difficult  or  impossible  to  integrate  over  the  intervals  of 
functions  where  the  slope  approaches  infinity.  In  the  current  version  it  is  still  diffi¬ 
cult  but  if  a  large  number  of  iterations  are  used  and  the  precision  parameter  is  large 
enough,  an  accurate  answer  will  be  obtained.  Test  Case  1  (Appendix  C)  of  the  current 
report  illustrates  this.  The  slope  of  this  function  approaches  infinity  at  the  start  and 
end  of  this  interval.  By  allowing  11  iterations  and  holding  quantities  to  10  places,  an 
answer  accurate  to  5  places  was  computed. 

Breakdown  will  definitely  occur  when  the  integral  is  close  to  zero,  resulting 
from  an  integrand  taking  on  both  positive  and  negative  values  over  the  interval  of 
integration.  A  function  such  as  this  should  be  integrated  in  parts. 

Input  for  the  program  has  been  changed  slightly  to  make  this  division  of  the 
integration  interval  more  amenable  and  there  are  no  limits  on  the  number  of  parts 
used.  This  situation  reveals  itself  readily  by  the  erratic  behavior  of  DELTA  (Aj). 
Test  Case  2  (Appendix  C)  illustrates  what  happens  when  integration  is  attempted  over 
the  whole  interval  of  such  an  integrand  and  Test  Case  3  illustrates  the  proper  method 
of  splitting  the  integration  into  parts . 
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m.  INPUT  LIMITATIONS 


The  integrand,  f(x),  may  be  any  function  for  which  the  integral  exists;  also 
(-<*>)  <  a,b  <  (+«),  b^a. 

The  total  number  of  iterations  allowed  is  Jmax  (an  input  parameter)  where 
2  s  Jmax<  100.  Iteration  is  stopped  whenever  either 

A  or  J  =  Jmax 

While  Jmax  can  be  any  number  up  to  100,  it  should  not  be  too  high.  Each  itera¬ 
tion  roughly  doubles  the  machine  time.  Generally,  for  more  than  13  or  14  iterations, 
machine  time  is  measured  in  hours . 

The  integrand,  f(x),  is  provided  by  a  subprogram  of  the  FUNCTION  type  named 
FUN(X).  FUN  must  appear  on  the  left  of  an  arithmetic  statement  whose  right  is  the 
integrand.  The  subprogram  source  deck  must  be  preceded  at  compilation  time  by  a 
precision  parameter  card  identical  to  the  precision  parameter  card,  if  any  was  used, 
of  the  main  source  deck. 

IV.  INPUT  AND  OUTPUT  FORMATS 

Input  Card  Format 

The  following  quantities  are  specified  in  the  input  (floating  point  numbers  are 
in  lower  case,  fixed  point  numbers  in  upper  case). 

Card  1 

Columns  1-80  contain  alphanumeric  identification.  Tt  is  reproduced  in  the 
output  but  does  not  affect  the  processing. 
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Card  2 


Columns 

1-36 

Integrand  Specification  Card.  Contains  alphanumeric 

information  to  identify  integrand  in  output.  This  does 

not  affect  processing. 

Card  3 

Format 

Columns 

1-12 

a 

E  12.6 

Columns 

13-24 

b 

E  12.6 

Columns 

25-28 

A 

F  4.4 

Columns 

29-30 

^max 

I  2 

Following  card  3  may  be  any  number  of  cards  specifying  new  values  of  a,  b,  A, 
and  Jmay.  Input  format  is  the  same  as  card  3.  The  final  card  is  a  trailer  card  with 
zeros  punched  in  columns  29-30. 

Input  Deck  Assembly 

Input  decks  are  assembled  as  in  the  original  NUMIT-ONE  with  the  FUN(X)  sub¬ 
program  replacing  the  POINT  subprogram. 

Errors  in  Input 

The  input  error  check  is  the  same  as  before  except  that  in  the  amended  version 
Jmax  may  have  any  value  up  to  99.  A  check  is  made  to  see  if  Jmax  <  2. 

Output 

The  output  has  three  main  stages;  the  identification  is  punched  as  follows: 

(a)  Alpha  Identification  Card 

(b)  Program  name  (NUMIT-ONE) 

(c)  F(X)  =  Columns  1-36  of  Integrand  Specification  Card. 
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Next,  the  input  a,  b,  Jmax>  and  DELTA  1  are  punched  followed  by  the  results 
of  the  computation  with  the  following  data  given  for  each  iteration  (J  =  2,  ...  Jact): 

(a)  Iteration  number 

(b)  Value  of  integral 

(c)  DELTA 

The  final  value  for  DELTA  £  DELTA  1  or  Jact  =  Jmax  is  noted  by  the  term 
(PASS  N)  in  the  last  line  of  the  output,  where  N  =  1,2. . .  99  is  the  pass  number. 

If  more  cards  specifying  new  values  of  a,  b,  Jmax  and  DELTA  1  are  input,  the 
above  cycle  will  be  repeated  until  all  input  cards  have  been  processed.  The  final  part 
of  the  output  is  furnished  when  the  trailer  card  is  read.  The  number  of  passes  (N)  is 
punched,  followed  by  the  total  value  of  the  integral  for  all  passes.  The  last  card 
punched,  showing  satisfactory  completion  of  the  program,  is  "END"  and  the  program 
stops. 

A  set  of  output  listings  for  sample  problems  is  given  in  Appendix  C. 

V.  OPERATION 

Operation  is  the  same  as  before  except  a  trailer  card  is  now  used  to  stop  the 
computation.  No  switches  are  used. 
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APPENDIX  A 


Plow  Charts 


Qiim  _ Sum  +  QQWfN  " 

U 1 1 1  »._  *  u  1 1 1  x  f  L  \*-  / 

Store: 
f(I)  —  fun 


n 


1  r 


Point  No. 

Flow  Chart  No. 


Flow  Chart  A-l 
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*Note: 


Point  No. 

Flow  Chart  No. 


a  -IAN1-AN2J 

AJ  ”  I  am  I 


Count  — 
Count  +  1 


AN2  = 

(4  Sumfx  +  Sum) 


Flow  Chart  A-2 
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APPENDIX  B 


Source  Program 

In  the  program,  the  following  symbols  are  used  which  are  unexplained  in  the 
mathematical  description: 

N  =  counter  to  keep  track  of  number  of  passes 

DELTA  1  =  A  =  input  convergence  criterion 


AN1 

=  Fj_i  =  previous  approximation 

AN2 

=  Fj  =  current  approximation 

SUMFX 

2J_1 

=  £  f(a  +  [2m- 1]  Ax) 

m=l 

SUM 

=  f(a)  +  f(b)  +  2S  2  •  f  (a+2m  *  Ax) 
m=T 

COUNT 

=  counter  used  to  count  successive  times  Fj  =  0 
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*1010 _  _  __  _  _ 

C  - NUMF  R IC  T  NTF  GR  A  T I  nN*~lIS  I  WG  SlMPSuN-  RVcKWPlIL-  ^FRRr-TWMC 

nlMFNSIttN  F (3) ,P<3) *0(3) ,U(9) 

p('i)=i. 

P(2)=4. 

p(3)=ir 

0(1)=1. _ _  _ 

Q ( 3) r 1 • 
n=p. 

TnTAl.rO . 

C 

C  REAR  INPUT  AND  ALPHA  IHFNT IFIC ATTnN  CARO 

~REAni02  . . 

_ PllNCHlO? _ 

PUNCH103 

READ  1 14 »  ( U (  T  )  »I=1»9) _ _ _ _ _  _ _ 

PUNCH l! 5 •  Ui(T ) r I  =  l79) 

1  RFAnim  »  A  >R»nFLTAl» UMAX _  _ _ _ 

C 

C _ TEST  INPUT  FpR  VALIDITY t  IF  pKAY  PUNCH  NFXT  PART  cF  PllTPUT _ 

IF(JMAX-2)?Gr26r27 

26  IF  (N)?Sf?5r?9_  _ _  _  _ _ 

25  PUNCH106~ 

ST  aP 

*27  IF(R-A)2?'r2R“r2? 

2A  PUNCH104 _ __ _ 

STnP 

22  IF(HFLTA1 123*24*23 _  _  _  _ _ _  __ 

24  PUNCH! 05 

STnP  _ 

23  PUNCH 109 *  B  * JM A X 

_ * _ PUNCH  110  *  A  *OFLTA  t _ _ 

c 

C  INI  TI  ALI7ATlnN_FnR  START _  ___  _ _ 

CpUnTsO • 

SUMrO  • _ _  _  _ _ 

StJMFXrO. . 

NrN+1 _ _ _ __ _ 

j=i 

K?=2  _  _  .  _ _ _ 

21  OFLYr  <R-A  )/<>.** J) 

IF(  J-l>4*_2*4  ___ 

?  On  3T— 1V3 

A 1=1-1 _ 

X  =  A  +  A I *OFLX 
F(I)rFUN(X) 

SUM  =  SUM+G ( T ) *F ( I ) 

3  SUMFX=SUMFX+P  (  I  )  *F  ( I  ) 

c  first  approximation  nF  integral 

ANlrOFLX*SUMFX/3 . _ _ _ 

U=J+1 
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no  n  n  n  r> 


<-.n  lu  PI 

C  CALC  ML  ATfS  FlINCT  InN  MinWAY  RFTWFFN  CURRENT  PnlNTS 

4  mi MFX=n# 

On  *=  1=1 *K2  _  _  _ __ _ _ _ _ 

AT  =  7*1-1 
X  =  A  +  A  T  *nFL  Y 

5  SIJlviFX  =  SI  1N.FX  +FUM  (  X  ) 

r.  NEXT  APPRnXTMATInN 

AN2=nFLX*  (  4  *  *SUMFX+SUM  )/3* 

IF( ANP)f>* 12*8 

8  OFLTA=ARSF( (  AN1 -AN2)/AN?*H10.  ) 

TF ( HFLTA-DFl  TA1  )17*17#18 

18  IFM-JMAY) IP* 17*17 

19  AN l r A NP 

SlJtf  =SUM  +P.*SUMFX 
PUNCH 1 1 2  *  J  *  AN 1  *  HFLT  A 

C  TNITT  ALT7ATTiiN  FnR  NFXT  ITERATION  I  nnp 

J=J+  1 
K2=P*KP 

C  BACK  Tu  TTFPATItiN  TF  VALUE  HAS  NnT  OnNVFRGEO  FNnlJGH 

Gu  Tu  21 

PUNCH  FINAL  ANSWFR  FcR  THIS  PASS  ANH  RfTuRN  FtxR  NFXlHlNPTiT  CTARf) 
17  PUNCH  1  13  *  J  *  N  *  ANP  *  DELTA 

TnTAl  -TnTAl  +AN2  . v' 

(;n  Tu  1 

FnR  HANOI  INC  7FRu  AMSWFRS 
17  CnU»vT=CmiNT+1  . 

T  F  (  CnUNT-7  •  ) 1 8  *  17  *  1 7 


PUNCH  THF  TnTAL  ANSWER  AND  HALT 
29  PUNCH  1 18*N*TnTAL. 

STiiP  _  _  _ _ _ _ _ 

lu  l  FORMAT  (  7F  12  #8*F4.4*I7) 

107  FnRMAT(anH 

1  ) 

103  FnH^AT(/l5X*4aHNUMIT  nMF  (uNF-HIM.  INTFGR  AT I nN *  SlMPStfMS  RULF)  ) 
11,4  FtiRN.-AT  (  1.  S  X  *  1 7HC HFCK  INPUT  -  A=R) 

105  FWRV  A  T  (  1 5X  *  7  2  HC  HFCK  INPUT  -  OFLTAIrO) 

108  FnR*.  aT(  15X#3PHCHFCK  INPUT  -  JMAX  miTSinE  FLANGE) 

1 OP  FaRi*AT ( 35X  *  1 OH INPUT  OATA/POX  *  PHAzFj  2,8* 11 X*7HJ(NAX)  =  * IP) 

110  FnRr  AT(7  0X*PHA  —  *F1P«8*  1  IX*  7HDFL.T  A1=*E10*4//1  5x  *  ToH  I  TFR  ft  T I  nN  *  1 1 X  * 

1 8HINTFGPAL *  1 3X  *  SHOFLT A ) 

IIP  F nH  v AT (IPX* T P *  1 PX *F 17 . 8* 1  OX  * F8 . 3 ) 

J 13  FnR  aT(19X*T7*SH(PASS*T3*1H) * 3X * F 1 P . 8  *  1  OX  * F8 . 3//// ) 
j.14  FnR  AT(PAS) 

11S  FnR  AT(1SX*7HF(X)  =  *9a5*//) 

118  FuR* AT ( 3PX *  1 7HNn .  PASSFS  —  I  3/P9X  *  1  8HTuTAl  INTEGRAL  =  Ft  P*5/39YT 
1 3HFNO ) 

FNO 
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APPENDIX  C 


Test  Cases 

Running  time  is  still  highly  variable,  but  in  general  the  time  is  now  30  to  50  per¬ 
cent  faster  than  in  the  old  version  as  will  be  noted  in  the  test  cases.  Also,  accuracy 
has  been  increased  since  quantities  are  now  held  to  10  significant  figures. 

Note  the  integrands  of  Cases  3  (a)  and  3  (b)  are  identical  to  Case  2  and  taken  to¬ 
gether  duplicate  Case  2.  However,  Case  3  produced  a  more  accurate  answer  in  a 
small  fraction  of  the  time  for  Case  2  because  the  Case  2  integral  took  on  positive  and 
negative  values  during  the  iteration  attempting  to  approximate  a  zero  answer  while 
Case  3  integrated  nonzero  portions  and  summed  them. 

Running  time,  as  in  the  original  report,  is  still  dependent  on  the  formula 

t  =  K  •  2 Jact 

Table  I  gives  the  new  running  times  for  the  test  cases. 


TABLE  I 


Test 

# 

Integrand 

a 

b 

3act 

t 

t  new 

.  ,  ,  x  100 

t  old 

K  (min) 

1 

V  (1-x3 )  (2-x) 

-i 

1 

11 

12'50" 

51.5 

6.27xl0-3 

2 

sin(x) 

0 

2n 

11 

16'50" 

64.9 

-3 

8.22x10 

3  (a) 

sin(x) 

0 

n 

5 

1* 

- 

8.22xl0-3 

3(b) 

sin(x) 

n 

2n 

5 

1' 

- 

8. 22x1 0-3 
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TEST  CASE  1 


INPUT  DATA 

SOURCE  STATEMENT  OF  SUBROUTINE 


*1010 

FUNCTION  FUN ( X ) 

FUN  =S0RTF( (l.-X*X)*(2.-X) ) 
RETURN 

END 


OUTPUT  DATA 


TFST  FUNCTION  1 


NUMIT  ONE  5KONE-OIM.  INTEGRATION f  SIMPSONS  RULEa 
F%Xn  tt  SQRT%%l.-X*Xn%2.-Xnn 


INPUT  DATA 

B# 

. 100000F&01 

J%MAXn#13 

Att- 

. 100000FA01 

DELTaT1«  . 

1TTO0E-02 

ITERATION 

INTEGRAL 

DELTA 

2 

.209138F&01 

9.838 

3 

.21 6390F&01 

3.351 

4 

• 21 8943F&01 

1.168 

5 

.219843FS01 

.409 

6 

•220161F&01 

.144 

7 

.220273FA01 

.050 

8 

.220312FS01 

.017 

9 

.220326FA01 

.008 

10 

. 220331F&01 

.002 

11%PASS  in 

.220333F&01 

0.000 

_  NO.  PASSES  tt  1 _ 

TOTAL  INTEGRAL  «  2.20333E-00 
END 
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TEST  CASE  2 


INPUT  DATA 

SOURCE  STATEMENT  OF  SUBROUTINE 


*1010 _ 

FUNCTUJN  FUN(X) 

_ FUN  r  SINF(X) 

RETURN 

END 


OUTPUT  DATA 


TFST  FUNCTION  2 

NUMIT  ONF  %ONE— DIM •  INTEGRATION*  SIMPSONS  RULED 

F%X«  #  SlN%Xn  


INPUT 'DATA 

B U  .628319F&01  J95MAXnttll 

A#  ,000000F-'59  OELTaI#  TTOOTjE-TT? 


ITERATION 

INTEGRAL 

delta 

3 

-.628319F-09 

100.000 

4 

-.523599F-09 

P.0.000 

5 

■ 641408F— 09 

181.632 

6 

- • 615229F— 09 

204.255 

7 

-.187841F-08 

67.247 

8 

— • 144971F— 08 

29.571 

9 

• 397608F— 09 

464.609 

10 

-.443300F-08 

108.969 

11SKPASS  in 

-.274807F-08 

61.312 

NO.  PASSES  tt  1 
TOTAL_  INTEGRAL  tt-2 . 74807F-09 

END 
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TEST  CASES  3(a)  and  3(b) 


OUTPUT  DATA 


TEST  FUNCTIONS  3(a)  and  3(b) 

NUMIT  ONF  KONE-DIM.  INTEGRATION*  SIMPSONS  RULEn 

F%XH  #  SlNSiXn 


Btt  .314159F801 

INPUT  DATA 

J55MAXH411 

Att  . 000000F—99 

DELTA 14  . 

1000E-02 

ITERATION 

INTEGRAL 

“DELTA 

2 

.200455F801 

4.481 

3 

. 200026F&01 

.214 

4 

.200001F&01 

.012 

5SKPASS  in 

•200000F&01 

0.000 

INPUT  RATA 

R«  .628319F801  J%MAXn#ll 

A#  . 314159E801  DELTA14  .1000E-02 


ITERATION  INTEGRAL  DELTA 

2  -.200455F&01  4.481 

3  -.200026F&01  .214 

4  -.200001F801  .012 

5%PASS  2n  -• 200000F801  0.000 


NO.  PASSES  n 

2 

TOTAL  INTEGRAL  U 

O.OOOOOF-99 

END 
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APPENDIX  D 


Glossary 

Fj(a,b) 

AXj 

f(x) 

CJ 

DJ 

J 

Aj 

A 

a 

b 

m 

^max 

Jact 

K 


the  approximation  of  the  integral  for  the  Jth  iteration  over  the 
interval  (a,b) 

the  mesh  interval  for  the  Jth  iteration 
the  integrand 

the  sum  of  the  values  of  the  integrand  over  the  interval  (a,b) 
whose  coefficients  are  2 

the  sum  of  the  values  of  the  integrand  over  the  interval  (a,b) 
whose  coefficients  are  4 

used  to  denote  iteration  number 

convergence  number  calculated  for  the  Jth  iteration 

input  convergence  criterion  number 

beginning  point  of  the  integration  interval 

end  point  of  the  integration  interval 

used  as  an  index  in  summation 

maximum  number  of  iterations  permitted 

actual  number  of  iterations 

constant  used  in  estimating  running  time 
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